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Abstract. We propose new semi-implicit numerical methods for the integration of 
the stochastic Landau-Lifshitz equation with built-in angular momentum conservation. 
f-H ' The performance of the proposed integrators is tested on the ID Heisenberg chain. 

For this system, our schemes show better stability properties and allow us to use 
\ considerably larger time steps than standard explicit methods. At the same time, 

these semi-implicit schemes are also of comparable accuracy to and computationally 
much cheaper than the standard midpoint implicit method. The results are of key 
importance for atomistic spin dynamics simulations and the study of spin dynamics 
beyond the macro spin approximation. 
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1. Introduction 

Dynamics of magnetic materials have been theoretically studied for many years starting 
from the seminal work by Landau and Lifshitz [1] (see, e.g., the monographs [2j [31 H]). 
The current interest in this area is rapidly growing due to new important fields of 
applications such as spintronics [5] and laser-induced ultrafast spin dynamics [61 [7J El [9] . 
In many situations such as the interaction of domain walls with pinning centers [ID], 
there are atomic-scale inhomogeneities which require multiscale simulations bridging 
macroscopic and microscopic lengths [TT]. In [121 03] a method of ab initio spin 
dynamics was suggested relating first-principle electronic structure calculations with 
Landau-Lifshitz-type dynamics of classical spins within the framework of the rigid-spin 
approximation. 

Thus, Atomistic Spin Dynamics (ASD) simulations are important from many points 
of view. To do calculations at finite temperatures, there are two main approaches: a 
generalized Nose-Hoover (Bulgac-Kuznecov) thermal bath or the Langevin (stochastic) 
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dynamics [13]. The first method has fictitious dynamics, and hence it can be used to 
simulate equilibrium properties only. The Langevin spin dynamics with first-principle 
magnetic interaction parameters has recently been implemented [H] and applied for 
simulating dilute magnetic semiconductors [T5] and spin glasses [IB] . Langevin spin 
dynamics are also used as a phenomenological simulation tool, not connected with first- 
principle theory. An implementation of this type was reported in [T7j and applied to 
laser-induced magnetization dynamics [T8] . 

The heart of Langevin spin dynamics simulations is integration of the stochastic 
Landau-Lifshitz (SLL) equation for each atomic spin. This equation is non-linear and 
analytical solutions for interacting systems exist for two spins only. In systems of interest 
for applications the number n of spins is typically of order 10 6 and the integration should 
be done numerically. Due to the interactions, one has to solve a system of 3n coupled 
non-linear equations. To compute quantities in equilibrium, this very large system 
should be simulated over long time intervals, usually from 10 fs to 1 ns. This is a 
challenging computational task. 

Thus, ASD requires effective numerical integrators for the SLL equation. Due to 
the large system size and long simulation time, such numerical methods should be, on 
the one hand, sufficiently stable and on the other hand very fast. The latter rules out 
the use of fully implicit integrators such as the implicit midpoint (IMP) scheme (see 
its application for Langevin spin dynamics, e.g., in [19]). Despite its superior stability 
properties which allows large step sizes, typically 10 fs, IMP is slow in practice since 
the implicitness requires solution of 3n non-linear coupled equations at every time step. 
Langevin spin dynamics simulations have often been based on the Heun method [HI [17] , 
which has the advantage of being fast in terms of the number of operations per time 
step. However, this method has poor stability properties requiring a relatively small 
step size, typically ranging from 0.01 fs to 1 fs, depending on the implementation. We 
also note that since the accuracy of the first-principle magnetic interaction parameters 
is limited to 10%, the accuracy of numerical methods is, to some extent, less important 
here than their stability (in the sense of the ability to use larger step sizes for long time 
simulations). Hence, both the standard implicit and explicit numerical integrators are 
not optimal for ASD and it is desirable to develop a numerical method that is both 
stable and fast. Also, ASD simulations are often used to study systems with different 
interactions and/or different symmetries. Therefore, in addition we should require from 
numerical integrators for ASD to be universal in their implementation. Such a method 
is proposed in this paper. 

As is known from the deterministic ([20j ETJ [22] and the references therein) and 
stochastic [231 I2H [25] numerical approaches, to numerically integrate dynamical systems 
over long time intervals with relatively large step sizes, it is advisable to preserve 
geometrical properties of the continuous dynamics. Therefore, one should construct 
and use geometrical integrators for ASD. 

In the case of the deterministic Landau-Lifshitz (LL) equation, there are geometric 
integrators [211 1221 [261 EZ] that are both stable and fast. Usually, these schemes are 
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semi-implicit. Unlike IMP, a semi- implicit method requires only the solution of 3 linear 
coupled equations for each spin individually. However, implementation of these methods 
depends on symmetry and interactions in a system under consideration, which makes it 
difficult to use them for models with arbitrary lattice structures. 

Further, semi-implicit methods for the deterministic LL equation are also considered 
in the review [28]. Being based on IMP, they have the potential to combine stability 
and low computational costs like the geometric integrators but with the advantage of a 
universal implementation. In this paper we use the idea of semi-implicitness to derive 
new numerical methods for Langevin spin dynamics simulations, which are both stable 
and fast and allow universal implementation. In particular, we show that, due to the 
enhanced stability, our semi-implicit integrator (named SIB) allows time steps by a 
factor of 10 -j- 10 3 larger than the standard Heun method. 

This paper is organized as follows. In Section |2} we formulate the problem in 
mathematical terms, introduce the necessary notation and examine the conservation 
properties of the SLL equation. In Section [31 we propose two new semi-implicit 
methods (SIA and SIB) and recall the Heun scheme and IMP. Both SIA and SIB 
intrinsically preserve the length of individual spins while SIB (like IMP) also possesses 
other conservation properties in the deterministic case. The later is apparently the 
reason for the superiority of SIB which is the numerical method of our choice for ASD. In 
Section HI we present some results of numerical experiments. We first test the considered 
numerical methods in the deterministic case without damping, using a simple system 
of two interacting spins. Then the ID Heisenberg chain is used as a test system for the 
stochastic case. In the last section, we draw conclusions and recommendations for future 
work. Two appendices are included to provide some auxiliary knowledge of stochastic 
numerics and about ergodicity of the SLL equation. 

2. Mathematical model 

In this section we formulate the problem in mathematical terms and introduce the 
necessary notation. In addition, we discuss why we use the Stratonovich interpretation 
for the stochastic LL equation. Finally, we examine the properties of the solution of the 
equations under study. 

The (deterministic) Landau-Lifshitz equation in dimensionless variables can be 
written in the form: 

__ = -X i x B l (X.) - aX i x \X l x B*(X)] , z = 1, . . . , n, (1) 

where n is the number of spins, X 1 = (X*,X*,X*) T are three-dimensional column- 
vectors representing unit spin vector^] and X = (X lT , . . . ,X nT ) T is a 3n-dimensional 
column- vector formed by the X 1 ; B l is the effective field acting on spin i; a > is 

X In the paper, we follow the standard notation of the theory of stochastic differential equations and 
use capital letters to denote solutions of differential equations while we use small letters for the initial 
data and for corresponding "dummy" variables. 
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the damping parameter. In (TTJ) the time is normalized by the precession frequency 
ujg = 7-B, where B is some reference magnetic field strength, and the effective field 
B = (B lT , . . . , B nT ) T is also normalized by B and is given by 

£(x) = -Vtf(x), (2) 

where H is the Hamiltonian of the problem. Then 

B*(x) = -Viff(x), 

where Vj is the gradient with respect to the Cartesian components of the effective 
magnetic field acting on spin i. 

For atomistic spin dynamics, the most important contributions to the Hamiltonian 
are the Heisenberg exchange for the interaction between the spins H ex , the Zeeman 
energy for the interaction with an external field H ext , and the uniaxial anisotropy H ani 
defining a preferential direction of the spins. Therefore we consider here the following 
Hamiltonian for our problem: 

H = H cx + i? ext + -Hani, (3) 



where 



ifex(x) = JjjXW, ffext(x) = -Bp ^ 

i^j i 
\2 



X 



F ani (x) = K^i 



x e K 



Here Jy are the exchange parameters, 5o is the uniform external field, K is the strength 
of the anisotropy, and ck is a unit vector that defines the anisotropy axis. Note that with 
these contributions to the Hamiltonian the effective fields B l are linear in x. In realistic 
materials usually |Jy| ^> \B \ 3> \K\. For the exchange parameters themselves, typically 
J%a+i) ^ Ji(i+j), J 1 > 1) i- e -; a U spins interact with each other but the nearest-neighbor 
interactions dominate. Since all the spins interact, Eq. (pQ) involves simultaneous solution 
of a 3n system of non-linear equations. Due to the interactions between the spins, each 
effective field B l is time-dependent and Eq. (JTJ) has in general no analytical solution. 
As a result, efficient numerical methods are required to study spin systems. In turn, 
the time-dependence of the effective field is usually considered as the main source of 
instability in the numerical integration. 

In order to perform spin dynamics at finite temperature, fluctuations are included 
according to the Brownian motion approach for spins by adding fluctuating torques to 
Eq. (CD) [291 130]. The stochastic Landau-Lifshitz (SLL) equation is then given by 

dX i 

— = -X l x (B l (X) + V) - aX l x [X 1 x (B l (X) + b 1 )} , (4) 
i = 1,.. .,n, 

where the fluctuating magnetic fields b l are uncorrelated Gaussian white noises 
interpreted in the sense of Stratonovich and 

(&{(*)) = 0, (b\(t)bi(0)) = 2D6 ij 6 lk S(t), z = l,...,n, (5) 



Stable and fast semi-implicit integration of the stochastic Landau- Lifshitz equation 5 



with (•) denoting ensemble averages and I, k = x,y, z labeling the Cartesian coordinates 
while D is the strength of the fluctuations. According to the fluctuation dissipation 
theorem, we choose 



where X is the (non-normalized) magnetization of each spin. 

Note that (Jl} is a differential equation with multiplicative noise which requires from 
us to specify in which sense we interpret the stochastic equation |31j. As said above, we 
use here the Stratonovich interpretation following [22]. This choice can be motivated as 
follows. First of all, the Stratonovich interpretation (contrary to any other one and, in 
particular, to the Ito interpretation) leads to preservation of the individual spin length 
(see (11 01) below) by (J4]), which is very important to model spin systems (see also a similar 
discussion in pj5]). Further, it is natural to model a perturbation of the Landau-Lifshitz 
dynamics by Gaussian noise with a finite bandwidth spectrum (i.e., by a colored noise 
[31J), possibly with a very short correlation time. The white noise b(t) in (jl]) has zero 
correlation radius (see ([5])) and a spectrum with infinite bandwidth. This noise is a 
convenient idealization which can be viewed as an approximation of the colored noise 
with short correlation time. Indeed, if we consider a sequence of solutions X n (t) of the 
equations X l n = -X l n x (^(XJ + tf n ) - ai; x [X l n x (5 i (X n ) + 6* )], where b n {t) is a 
sequence of Gaussian processes which correlation functions that go to the 5-function as 
n — > oo, then X n tends to the solution X of (jl]) if it is interpreted in the Stratonovich 
sense [32J Chapter 2], [33] Chapter 5]). We also note in passing that one can model a 
Gaussian colored noise by the Ornstein-Uhlenbeck process [31] which can be substituted 
in (jl]) instead of the white noise b(t). It could be of interest to study the influence of 
the correlation radius on the stochastic Landau-Lifshitz dynamics. We do not pursue 
such questions in this paper but remark that effective numerical methods for differential 
equations with colored noise are available in [HJ |32] which can be adapted to the SLL 
equation with colored noise. 

Since we will exploit some results from stochastic numerics [24J which in turn follows 
the standard theory of stochastic differential equations, it is convenient to re-write the 
SLL equation (Jl]) in differential form [31]: 



where W\t) = (W*(t), W*(t), W*(t)) T , i = 1, . . . , n; W*(t), WMt), W l z (t) } i = l,...,n, 



are independent standard Wiener processes; a,(x), x G M. 3n , are three-dimensional 
column-vectors defined by 



dX l = X i x Oi(X)dt + X i x a(X' 1 ) o dW\t), 
X l (0) = Xq, \x l \ = 1, i = 1, . . . ,n, 



(7) 



Oi(x) = -£ l (x) - ax 1 x £ l (x) ; 
and cr(x), x G M 3 , is a 3 x 3-matrix such that 



(8) 




(9) 
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for any y G M 3 . Note that the symbol 'o' in Eq. ([7]) means that the corresponding 
stochastic integral is interpreted in the Stratonovich sense [31]. We recall [32] (see also 
[3T|l33]) that the Stratonovich stochastic integral can be defined as the mean-square limit 
of the middle Riemann sums, which, in particular, makes it evident why the midpoint 
scheme (see (1151) below) satisfies the Stratonovich calculus. 

Let us consider some properties of the solution to ©-([9]). First, the length of each 
individual spin is a constant of motion, i.e., 

\X\t)\ = 1, i= l,...,n, t>0. (10) 

Indeed, we have 

d-\X l \ 2 = X i dX i 
2 1 1 

= X i [X* x o<(X)] dt + X 1 [X { x a(X*) o dW^t)] = 0. 

Other general conservation laws of ([7j)- (Q and also of ([I]) do not exist. However 
when we restrict ourselves to realistic systems, we have the damping coefficient a C 1. 
This means that, in practice, solutions of ©-(JHJ) are, in a sense, close to the deterministic 
solutions of (CQ) with a = 0. Hence the precessional motion can usually be considered 
as dominant. In turn, the largest contribution to the precessional motion is due to 
the exchange interaction. Therefore, it is relevant to examine the conservation laws for 
a = 0. Since the Hamiltonian has no explicit time-dependence, energy is conserved for 
this case. Further, when only Heisenberg exchange is included we have for the total 
spin: 

dX l 

J2^r = Y, J v xl x x ° = H J « ( x * x xj+xj * = o (ii) 

i i^j i>j 

since Jjj = Jj^. We recall that the orientation of individual spins is time dependent, 
which makes the effective field acting on each spin time dependent due to the exchange 
interaction. However, at the same time, the symmetry of the exchange interaction 
ensures that the total spin is time- independent. Therefore the conservation of total spin 
is an important property for stable numerical integration of the exchange interaction. 
By the same arguments, when an external field is added, the total spin will precess in 
the external field: 

= < 12 > 

i i 

For this case, the length of the total spin is a constant of motion, as well as the 
component of the total spin along B . Hence the energy is also conserved but the 
transversal components of the total spin with respect to B Q oscillate in time. When 
anisotropy is included, there are no conservation properties associated with the total 
spin. Finally, ergodicity of the solution to is a relevant property. This is discussed 

in Appendix A. 
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3. Numerical methods 

In this section we consider numerical integrators for the stochastic Landau-Lifshitz 
equation (J7])- ([9]). We first recall two existing numerical methods, one of which is 
explicit (the projected Heun scheme) and the other implicit (the midpoint scheme). 
Both are unsatisfactory since either they violate conservation laws (HeunP) or they are 
computationally very expensive (IMP). Therefore, in the main part of this section we 
present the two newly developed numerical methods (SIA and SIB). These methods are 
called semi-implicit and aim at combining the advantages of the existing explicit and 
implicit schemes. 

As it is known from the deterministic ([201 [211 122] and the references therein) 
and stochastic ([231 I2H [25]) numerical approaches, to achieve accuracy in long-time 
simulations (e.g., for computing ergodic limits) it is advisable to preserve the structural 
properties of the continuous dynamics by the approximating discrete ones. Then it is 
important to consider not only orders of convergence but also structural properties of 
numerical integrators for the SSL equation. Both convergence and structural properties 
of the schemes presented are discussed in Section 13.31 

Throughout we use (for simplicity) a uniform discretization of a time interval [0, £*] 
with step size h = t+/N. The value at the initial step is X l = x l Q , i = 1, . . . , n, and XI, 
i = 1, . . . , n, denotes the approximate solution X l (t k ), i = 1, . . . , n, to the SLL equation 
at time t k , k = 1, . . . , N. 

3.1. Existing explicit and implicit numerical methods 

3.1.1. "Heun + projection (HeunP)". The Heun method can be seen as a predictor- 
corrector method. Its prediction step, which we denote by X k , is the Euler 
approximation. The standard Heun method should be adjusted by an additional 
projection step which is needed to ensure that the length of each individual spin remains 
constant. For the SLL equation (EJ)-^nj), the HeunP method reads 

X* =Xl + hXt x a t (X k ) + h^Xl x o{Xi)& +x , (13) 
i = l,...,n, 

X*Ui = X l + \ [ X l x + $ x «*(**)] 

+ ^ [XI x a{Xl)C k+1 + XI x a{Xl)C k+1 ] , 

X k+ i = X k+1 /\X k+1 \, i = 1, . . . , n, 
k = l,...,N, 

where X k = {Xf, . . . , XfY~ ^ = #J X , 5 3 = h 2, 3, i = 1, . . . , n, 

k = 1,...,N, are independent identically distributed (i.i.d.) random variables which 
can be distributed, e.g., as 

P{& = ±1) = 1/2 (14) 
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or ^ J ~ A/"(0, 1). This indicates that the are i.i.d. Gaussian random variables with 
zero mean and unit variance. In Eqs. f lT3|) we explicitly added i — 1, . . . , n to emphasize 
that first Xk has to be calculated for all spins, before X^+i is computed. We come back 
to this point in the numerical experiments (Section @J . 

3.1.2. "Implicit Midpoint (IMP)". Contrary to the HeunP method, IMP (see, e.g. [2S 
p. 45])) is implicit. For the SLL equation (EJ)- (EJ) , IMP reads: 

Xi +1 = Xl + h 3±3±l x ai (15) 
k = l,...,N, 

where = (C k + v Ck+v &V ! J = 1, 2, 3, i = 1, . . . , n, fc = 1, . . . , N, are i.i.d. 
random variables which can be distributed according to, e.g., ffl4|) . Alternatively, we 
can choose Ck being distributed as the £/, defined below (see [23, 1211 )• C ~ -^(0, 1) 
be a Gaussian random variable with zero mean and unit variance. We define 

C, ICI < A h , 

U { A h ,C> A h , (16) 
-A h , C < -A h , 



where A h = y/2\ kih\. We note that if one takes Ck ~ 1), IMP can, in general, 
diverge (see a counter-example in [23| l2lj). 

5.^. New semi-implicit numerical methods 

Here we propose two new semi-implicit integration schemes, simply called semi-implicit 
A (SIA) and semi-implict B (SIB). In the spirit of the review [28] , they are called semi- 
implicit since they require only to solve n or, In in the case of the SIB scheme, linear 
3x3 systems at each time-step, which can be done analytically. The starting point 
for derivation of the semi-implicit methods is the IMP scheme. To reduce the degree of 
implicitness, we replace Xfc +1 in the argument of aj and a in IMP by a predictor X k . As 
a consequence, resolving the implicitness at each time step is simplified (in comparison 
to IMP) to solving a linear 3x3 system per spin that is independent of the interactions 
between the spins. The difference between SIA and SIB is the choice for X k . Both 
semi-implicit methods have effectively the same computational cost as explicit schemes. 

3.2.1. u Semi-implicit scheme A (SIA)". Similar to the HeunP method, for the SIA 
scheme we take the Euler approximation for the predictor X k . The SIA method for the 
SLL equation reads 

XI = XI + hX\ x a,,(X fc ) + h^Xl x a(Xl)a +1 , (17) 
% = l,...,n, 
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k=l,...,N, 

where = (Ck+n Ck+n €k+i) > Cl'*' are i-i-d- random variables as in IMP (|T5|) (the 



same two possibilities). 

3.2.2. u Semi- implicit scheme B (SIB)". SIA can be viewed as a second iteration for 
the implicit equation due to IMP. As zero approximation of Xfc +1 , we took and then 
the second iteration was constructed so that the length of individual spins is preserved. 
One can see that the first iteration (or in other words the prediction step) of SIA does not 
preserve the spin length. We are therefore proposing the SIB method which keeps the 
spin-length conserving IMP structure at both iterations and, according to our numerical 
tests (see Section 0J), this modification is crucial for the performance of the semi- implicit 
schemes. 

The SIB method for the SLL equation reads 



k = l,...,N, 

where = {Ck+v Ck+v Ck+i) > are i-i-d- random variables as in IMP (115]) (the 
same two possibilities). 

Remark 1 One can continue the process and make several iterations for the implicit 
equation due to IMP, e.g., in our tests about 10 iterations were sufficient to resolve 
the implicitness up to the machine accuracy. However, in practice the use of several 
iterations would be too computationally expensive while SIB already demonstrates 
stability and accuracy comparable with IMP. 

3. 3. Properties of the methods 

We start by examining convergence of the methods presented in this section and then 
discuss some conservation properties. For completeness, in Appendix B we recall some 
generic facts about stochastic numerics [21]. 

All four methods considered in this section are of weak order one for both choices 
of the distributions of £jf (discrete and continuous). If C( 3 ~ A/"(0, 1), then HeunP is 
also of mean-square order 1/2. IMP, SIA, and SIB are of mean-square order 1/2 if £jf 
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have the cut-off Gaussian distribution (TTB]) . These convergence properties are proved 
using the standard results [2H Chapters 1 and 2]. In the deterministic case (i.e., D — 0) 
all four methods are of order two. 

Note that in this paper we limit ourselves to methods of weak order 1 and of mean- 
square order 1/2. The system ©-(El) has noncommutative noise (see the definition 
in, e.g. [2U p. 28]). Then mean-square methods of orders higher than 1/2 require 
simulation of multiple Ito integrals which is computationally expensive. It is possible 
to construct higher order weak methods for (ED)-© but, due to the multiplicative, 
noncommutative nature of the noise, they would be too complicated and they are 
not considered here. We also note that the problem with multiplicative noise can be 
circumvented by rewriting the SLL equation in spherical coordinates, for which the 
system is Hamiltonian and the noise becomes additive, but then numerical difficulties 
arise when the polar angle is close to or 7r. 

When a is small, the SLL equation (E])-© is a system with small multiplicative 
noise. In this case the weak-sense errors of all the methods considered in this section 
are of order 0(h 2 + a 2 h) |35j,[24, Chapter 3]. The smallness of noise can be further 
exploited to construct high accuracy but low order efficient methods following the recipe 
from [351 US- 

We now discuss conservation properties of the schemes. The HeunP method (Tl3|) 
has only one conservation property - norm-preservation which is due to the projection 
step. Heun without the projection step would conserve the total spin but then violates 
norm-preservation. Omitting the projection step also gives very poor results for the 
interaction with an external magnetic field. In practice the projection step can be 
exploited for error control. Energy is not conserved by HeunP when a = 0. HeunP has 
the advantage of being very flexible, its implementation is independent of the symmetry 
of the system and types of interactions used. The method is also fast since integration 
can be done for each spin separately. 

Due to the structure of IMP, the difference X l k+1 — X k is always perpendicular to 
XI + Xt +V Therefore (X* k + X* +1 )(X* +1 - X{) = and hence |X£ +1 | 2 = |X|| 2 , i.e., 
the length of each spin is exactly preserved by IMP without any need of projection. 
In the deterministic case with a = and under only the Heisenberg exchange, IMP 
conserves the total spin. The proof follows directly from Eq. (TIT]) with replacing dX l / dt 
by X\ +1 — X\ and X 1 by (X l k + X k+l )/2. The total energy conservation for the case 
of a = can be proven similarly. Preservation of all the main structural properties 
of the SLL equation by IMP comes at a cost. Since all spins are coupled, a system of 
3n non-linear algebraic equations has to be solved at each time step. This is a major 
limitation for application of IMP to atomistic spin dynamics, where the number of spins 
is typically of order n = 10 6 . Some further remarks on conservation properties of both 
HeunP and IMP in the deterministic case are given in [2TJ . 

The SIA method is very close to the HeunP method. However, unlike the HeunP 
method SIA preserves the constraint |X l (t)| = 1 exactly, without the need of projection. 
This follows directly from the observation that the norm conservation of each spin is 
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independent of the point at which a, and a are evaluated. Let us now look at SIA in the 
deterministic case with a = 0. Regarding total spin, the relevant symmetry property is: 

xi + xi +l x x[ + x[ x{ + x{ +l xi . xj i(] 

2 2 2 2 T ' 

which is violated since the Euler approximation for X k depends only on the orientation of 
the spins at the current time step (X^), but not on Xfc, whereas X l k+1 is also determined 
by the value X l k+1 itself. Owing to this difference, for a = the total spin cannot be 
preserved by SIA. Also, the energy is not a conserved quantity by SIA and the scheme 
introduces numerical damping. Hence SIA has the same conservation properties as 
HeunP, and it is of interest to investigate whether the built-in norm conservation is 
sufficient to improve stability properties. 

Unlike SIA, SIB has the norm-conserving midpoint structure for both X l k and X k . 
In the case of a two-spin deterministic system with a = we proved analytically that 
both energy and total spin are conserved quantities of SIB. Hence for this system SIB 
has the same conservation properties as IMP. At the same time, implementation-wise 
very little additional computational efforts are required by SIB compared to HeunP and 
SIA. Hence it is of interest to compare the performance of SIB with SIA, in particular 
to investigate the influence of preservation of norm-conservation and preservation of 
deterministic conservation laws on the stability properties of the methods. As our 
numerical experiments (see the next section) suggest, SIB outperforms SIA while SIA is 
only slightly better than HeunP. This observation implies, in particular, that the built-in 
norm conservation alone is not sufficient for obtaining superior numerical integrators for 
ASD and preservation of other structural properties of the SLL equation should guide 
one in constructing effective numerical methods. 

4. Numerical experiments 

In this section, we compare performance of the integrators introduced in the previous 
section using two model problems. In Sec. 14. 1[ we present some results of the experiments 
in the deterministic case without damping (i.e. a = 0), to illustrate the conservation 
properties of the numerical methods. In Sec. I4.2[ we consider the stochastic case using 
the ID Heisenberg chain as a test system. We show that the methods that preserve the 
deterministic integration laws give rise to a more stable integration for the stochastic 
spin dynamics. 

4-1. Two interacting spins 

In order to illustrate the conservation properties of the numerical schemes related to the 
deterministic precessional motion, we choose the simple case of two interacting spins 
with equal length \X 1 \ = \X 2 \ = 1. As a result of the exchange interaction, the spins 
rotate around a common axis, where the precession frequency is given by uj = 2 J cos 8/2 
with the angle 9 between the spins and the Heisenberg exchange parameter J. 
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First, we emphasize the relevance of simultaneously updating the effective field. Due 
to the interaction, the effective field acting on each spin is determined by the other spin. 
Therefore, when using a predictor-corrector method like HeunP, it is highly relevant 
to simultaneously update the effective fields after the prediction step before calculating 
the correction step. Hence, the correction step is computed taking into account that 
ai{X k ) depends on and not on X£ alone. Therefore, at each time step the effective 
field must be computed twice. By its design, a predictor-corrector method must be 
implemented in this way, otherwise it will, as a rule, become a scheme of lower order. 
Figure [1] shows the computed trajectory with and without simultaneous update for the 
HeunP method. To achieve a comparable accuracy without simultaneous update of the 
effective field, the step size should be decreased by a factor of 10 2 -j- 10 3 . 

In the four lowest panels of figure [1] we compare the considered integrators 
implemented with simultaneous update of the effective fields. For illustration purposes, 
a large step size is used (h = 1/16). For small times, all methods show reasonable 
agreement with the analytical solution, but IMP clearly has the best performance for 
this system. However, even IMP, which preserves the conservation laws instrinsically, 
introduces errors in the precession frequency. Since these errors do not effect the 
conservation properties of the methods, we do not consider them in detail. 

Next, we compare the conservation properties of the considered methods for the 
2-spin system. To this end, figure [2] shows the error in the total spin as a function of 
integration time. Both SIB and IMP exactly conserve the total spin, whereas HeunP 
and SIA have numerical dissipation. For clarity, only the ^-component of the total spin 
is plotted. The errors in the x and y-components of the total spin are much smaller 
since the numerical errors in the x, y motion of the individual spins cancel each other 
due to the symmetry. 

Despite the fact that SIA conserves the norm of each spin exactly, the numerical 
damping is slightly larger than for HeunP. Both their errors are strongly dependent on 
the initial condition. When the spins are almost parallel, HeunP has a larger numerical 
error than SIA since the projection step transforms a larger amount of transverse motion 
to longitudinal motion. In the case of figure [2] an initial condition with 8 = 120° is 
used, which is closer to anti-parallel motion and, therefore, HeunP has a smaller error 
than SIA. 

For this simple 2-spin system, the energy and total spin are directly related: 
(Xl + X 2 k f = (X l k f + (Xlf + 2X\X\ = 2 + E k /J. Hence both SIB and IMP conserve 
energy, whereas both HeunP and SIA dissipate energy. For larger systems with only 
nearest neighbor interactions, SIB conserves total spin and energy like IMP as well, while 
obviously SIB requires much lower computational efforts than IMP. The conservation 
properties of SIB can be proven analytically but this is beyond the scope of the present 
paper. 

In conclusion, the results of the numerical experiments with 2 interacting spins 
and a = show that both HeunP and SIA introduce numerical errors in the conserved 
quantities whereas SIB and IMP preserve the total spin and energy of the test system. 
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Figure 1. Comparison of the explicit HeunP, implicit IMP and semi- implicit methods 
SIA and SIB for the deterministic case a — 0. The trajectory of 2 interacting spins 
is shown by plotting the x components of the 2 spins and 1 z-component. Solid lines 
indicate the analytical solution. The upper panel shows that without simultaneous 
update of the effective field the integration is very unstable. IMP demonstrates the 
best performance. All methods introduce errors in the precession period tj — 2tt/J 
corresponding to initial condition. For the purpose of illustration, a large step size 
h = 1/16 is used. 



4-2. ID Heisenberg chain 

In this section we compare the semi-implicit integration schemes with the explicit and 
implicit methods in the stochastic case. The simplest model of classical interacting 
spins is the ID Heisenberg chain with nearest-neighbor interactions. For this system, 
an analytical expression for the mean energy per spin is available [361 13?] : 

This expression gives us a convenient way to check how accurately the temperature of 
the system is reproduced in simulations using the numerical methods from Sec. [3j Note 
that H — > — 1 + 1/n as the temperature T — > since we have normalized the energy 
with the number of spins n and the interaction energy of 2 spins 2JX l X 2 tends to 2J 



Stable and fast semi-implicit integration of the stochastic Landau- Lifshitz equation 14 
0.02 



-0.02 

c 

Q. 

W 

■3 -0.04 

-I— ' 

O 

H — 1 

-0.06 

O 
i_ 

LLJ 

-0.08 



-0.1 



-0.12 




time (units of t ) 
j 



Figure 2. Conservation of total spin for HeunP, IMP, SIA, and SIB. Shown is the 
error in the total spin for the same system as in figure [TJ Both IMP and SIB preserve 
the total spin up to machine precision, whereas SIA and HeunP introduce a numerical 
damping. Here tj = 2ir/J is the precession period. 



when the temperature goes to zero. 

The comparison of the HeunP method with the semi-implicit schemes for the 
temperature is shown in figure [3] for step size h = 1/32, damping a = 0.1, exchange 
parameter J = 1, spin length \X % \ = 1, and number of spins n = 100. The random 
variables used in the numerical schemes are simulated according to the cut-off Gaussian 
distribution f fT6|) . At a time step k the sample average Hk for the energy H per spin is 
computed as 

m=l 

where Xj^ are independent realizations of X*. obtained by a numerical scheme (see 
also Appendix B). The corresponding standard deviation au k is also computed. In the 
experiment an ensemble of M = 20 independent trajectories was used. The values 
plotted in figure [3j with the 95% confidence intervals determined by the standard 
deviation, were obtained after equilibrating the system for a time t a = 1024 tj, long 
enough for the system to be sufficiently close to equilibrium. Here tj = 2n/(2J) is 
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Figure 3. Temperature check of of the semi-implicit methods SIA and SIB compared 
with the explicit HeunP method. Shown is the mean energy per spin of the ID 
Heisenberg chain, as function of temperature, computed with the parameter values 
shown at the bottom. All the schemes demonstrate reasonable agreement with the 
analytical result (TUT 



the reference precession period for (almost) parallel spins. We find that both HeunP 
and the semi-implicit schemes show reasonable agreement with the analytical results, 
indicating that they obey the Stratonovich interpretation rule as expected. 

The next question is which method is more accurate. Figure |3] shows that SIB 
is consistent with the analytical solution at all data points. To the contrary, HeunP 
and SIA show slight discrepancies. To investigate this more accurately, we study the 
numerical error by varying the step size. For illustration, we used the lowest temperature 
k},T/(2J) = 0.1. The results are shown in figure HJ 

It is found that SIB outperforms both SIA and HeunP, and SIB remains stable 
down to only 4 steps per precession period. At such a large step size, SIA and HeunP 
are unstable though SIA performs slightly better than HeunP. Note that in physical 
units, with the exchange energy JX 2 = 1 mRy, X = 1/iBohr, 4 steps per precession 
period corresponds to a step size of about 20 fs. Hence, for SIB the step size is only 
limited by the precession period of the spins, and there is no need to decrease the step 



Stable and fast semi-implicit integration of the stochastic Landau- Lifshitz equation 16 



o 
co 

H — ■ 

& 

CD 

c 

CD 

CO 
CD 



O 
CD 



1.4 
1.2 
1 

0.8 
0.6 
0.4 
0.2 



-0.2 



- © - HeunP 
-a -SIA 

-a - SIB 



¥ 



i i 
I i 
i i 
i i 



i 



/ i 
/ i 

0**6 = 4 = ±_ 



- -e- 4 



0.05 



0.1 0.15 
step size h (units of t ) 



0.2 



0.25 



Figure 4. Stability of the semi-implicit methods SIA and SIB comared with the 
explicit HeunP method. Shown is the error in the mean energy as a function of the 
step size h for the lowest temperature considered in figure [21 fc&T/(2J) = 0.1. It is 
found that SIB remains stable up to 4 steps per precession period tj = 2ir/(2J), while 
SIA and HeunP become unstable and produce unreliable results. 



size to preserve the conservation laws accurately enough. This should be compared with 
the step size of 10 as which was reported in [13] , resulting in an enormous improvement 
of a factor 2 ■ 10 3 in the allowed step size. However, the mentioned implementation 
of ASD in [14] is based on HeunP without the simultaneous update of the effective 
field. As follows from figure [T] and figure HJ when the effective field is properly updated, 
HeunP also allows a larger step size. However, the increase is limited to about 2 fs 
for the system studied here. Compared to HeunP, SIA has only slightly better stability 
properties, which we attribute to the intrinsic norm conservation. The superior stability 
properties of SIB can apparently be explained by its built-in deterministic conservation 
properties. For the system studied here, SIB allows step sizes by about a factor of 10 
larger than HeunP and by about a factor of 5 larger than SIA. 

Let us now compare the performance of the semi-implicit methods with the full 
implicit IMP. The ID Heisenberg chain is not convenient for this purpose, unless we 
choose a very small number of spins. In addition, for this comparison stability is not 
the major issue since we already know that the step size of SIB is limited only by the 
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Figure 5. Comparison of the semi-implicit methods SIA and SIB with the full implicit 
IMP. Shown is the weak-order convergence of SIA, SIB and IMP schemes for the mean 
energy per spin. Both axes are logaritmic with base 2. For small enough step size, 
the slope gives the order of convergence. Surprisingly, both SIA and SIB are more 
accurate than IMP. Moreover, SIB shows a higher order convergence than IMP. Here 
tj = 2tt/(2J) indicates the reference precession period. 



precession period. Therefore, we are more interested in the intrinsic properties of the 
integrators that are independent of the system under study. Hence the relevant property 
here is the convergence of the semi-implicit and IMP schemes. To reduce computational 
costs of the experiment, we again use a system with only 2 spins. 

To experimentally observe the order of convergence, a small statistical error is 
needed. To this end, a combination of ensemble and time averaging was used. As 
before, for an ensemble with M trajectories, we let the system equilibrate for a time 
t a = 2048 tj. Subsequently, the equilibrated sample mean H k (see ( ED]) ) is calculated for 
a time i& = 6144 tj. The calculated values of Hk are then divided in P = 8 subsets of 
length L = tb/P = 768 tj and in each subset the time mean H p is computed. Eventually, 
the total mean H is the average of the time means over the P subsets and its statistical 
error A is estimated by two standard deviations of H p divided by yP, which gives half 
of the length of the 95% confidence intervals for H. The results are presented in figure [5j 
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Note that for this small system no instabilities appear in SIA, and this method 
shows the first weak-order convergence as expected. Surprisingly, SIB demonstrates 
a second-order convergence, which might be related to the fact that the energy is a 
conserved quantity for a = 0. This means that for the energy only numerical errors 
from the damping term show up, hence the convergence for the energy in the stochastic 
case might be better than the convergence for a general quantity. The small error for 
SIA at the one but smallest time step in figure |5] is caused by the change in sign of the 
error. The error values are given in Table HJ Here also the data for HeunP are provided. 
HeunP is not shown in figure |5] since it appears to be in the asymptotic regime only for 
the smallest time steps. We note that there is a sign change of the HeunP error, which 
is the reason for its small error at h — 1.564 x 10~ 2 . IMP is very costly for a large 
ensemble, therefore the two smallest step sizes were not computed. 



Table 1. The values of error in the mean energy e = H — H analytic and the 
corresponding statistical error A for the considered schemes. In each consecutive row 
the step size is smaller by a factor 2. 
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In general, the performance of SIB in the experiments has been better than SIA. 
Interestingly, despite the excellent stability of IMP, the accuracy of IMP in the stochastic 
case lags behind SIB and SIA. This is a good example of a situation when a method 
with better stability not necessarily has a better accuracy. It was also observed in the 
deterministic case with damping that SIB sometimes shows better accuracy than IMP. 
This implies that in the case of damped motion the numerical integration error of IMP 
can be larger than for SIB, as it is observed in the stochastic case. These results show 
that at least for the systems considered here, SIB has the same stability properties as 
IMP, but at considerable lower computational costs. 

In conclusion, we find that in the stochastic case the semi-implicit method B, with 
built-in deterministic conservation laws is more stable and has smaller numerical errors 
than both the SIA and the HeunP method. Surprisingly, in the stochastic case SIB is 
even better than IMP in terms of accuracy and convergence. SIA performs only slightly 
better than HeunP in the stochastic case, and from this we find that norm-conservation 
is not the most important criterion for stable numerical integration of the SLL equation. 
Hence, SIB combines the advantages of both HeunP and IMP, being both fast and stable 
as well as universal. For systems with only nearest neighbor interactions, SIB allows 
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step sizes by a factor of 10 larger than the popular HeunP scheme, and a factor of 
2 • 10 3 larger than the HeunP method without simultaneous update of the effective field. 
Since in practice nearest-neighbor interactions dominate, SIB is expected to be also 
advantageous for systems with more than nearest-neighbor interactions. 

5. Conclusions and Outlook 

In this paper we introduced two new semi-implicit integrators (SIA and SIB) for 
stochastic Atomistic Spin Dynamics (ASD) simulations. These schemes combine 
the advantages of the standard explicit projected Heun method (HeunP) and the 
fully implicit midpoint method (IMP). The semi-implicit methods are fast as explicit 
schemes since they require only the solution of 3 linear coupled equations for each spin 
individually and therefore they are effectively explicit. For stability, the most important 
conservation law is apparently the preservation of the total spin for the case without 
damping. Like IMP, SIB preserves this conservation law for the dominant interactions 
in the system and the stability properties of SIB are comparable with IMP. SIA, which 
has norm- conservation built-in but not the deterministic conservation laws, shows only 
slightly better stability than HeunP in the stochastic case. Therefore, we recommend 
the use of SIB for ASD simulations. 

Owing to the enhanced stability, larger step sizes can be used with SIB. From our 
numerical experiments we can conclude that the step size can be increased by a factor 
of about 10 compared to the explicit HeunP. For SIB, the step size is only limited by the 
precession frequency of the individual atomic spins in the exchange field, which allows 
for step sizes of about one fourth of the precession period which can be as large as 
20 fs. This value of the step size has to be compared with the 10 as that was reported 
for a standard implementation of ASD simulations [H], which is based on the HeunP 
method without the simultaneous update of the effective field. Hence, the factor 2 • 10 3 
improvement can be attributed to a proper update of the effective field and built-in 
conservation of the total spin for SIB. Interestingly, numerical experiments indicate 
that SIB can also be more accurate than IMP in the stochastic case. Further checks 
for the stochastic case, including larger systems, more complicated interactions, and 
correlations, will be discussed in a following paper. 

Future work should study the conservation properties of SIB in more detail in 
order to give a further explanation of its excellent behavior. It would also be of 
interest to obtain a method obeying conservation laws for systems with more complicated 
interactions (e.g. next-nearest neighbor, anisotropy). In addition, one might exploit the 
fact that the damping motion and the precessional motion are always perpendicular, 
which potentially can be used to design an integrator that exactly dissipates energy 
like in continuous dynamics. Another direction which we can pursue in future is to 
derive stochastic counterparts of the geometric integrators proposed in [21, 22J for 
deterministic Landau-Lifshitz equations. Though they lack flexibility to deal with 
models with arbitrary lattice structures, such geometric integrators are expected to 



Stable and fast semi-implicit integration of the stochastic Landau- Lifshitz equation 20 



be highly efficient when it is sufficient to include only nearest neighbor interaction in 
the stochastic model. Our method can also be of value for micromagnetic simulations 
and we expect that similar techniques can be exploited for other physical systems, where 
interactions between particles are governed by a global conservation law, e.g., systems 
based on diffusion equations such as the Schrodinger equation. 
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Appendix A. 

In this Appendix we discuss the ergodicity of the solution X(£) to ©-Q. For the 
solution X(i) of (0)- (EJ), we will also use the notation X x (t) to reflect the dependence 
on the initial condition X x (0) = x. Taking into account ([2]) and ([3]), we observe that 
the coefficients of ©-Q are smooth functions and due to ffTUl) they remain bounded 
for all t > 0. 

One can show [331 EH] that for D > and a > the process X(t) is ergodic, i.e., 
there exists a unique invariant measure /i of X and independently of x e M 3n there 
exists the limit 

lim (y>(X x (t))> = / <^(x) d//(x) = ^ (A.I) 



for any function tp(x) with polynomial growth at infinity. Indeed, the solution X(i) 
of (EJ)-© lives on the compact due to f lTUj) . Then to prove ergodicity, it is enough to 
show that there is sufficient mixing. When a = 0, the stochastic perturbation is only 
precessional and, in general (e.g., for constant B) the process X(t) is not ergodic. When 
a > 0, the stochastic perturbation acts in all the directions on the spheres \x % \ = 1 and 
so ensures a mixing sufficient for the ergodicity. 

We also recall the ergodic theorem, which gives the equivalence between the 
ensemble and time averaging: 

t 

1 



lim - / (p(X. x (s))ds = (p er9 almost surely, (A. 2) 

t->oo t 



o 



where the limit does not depend on x. 

Further, the invariant measure associated with the solution X(t) of ©-Q is 
Gibbsian with the density 

p(x) cx exp(-/3tf (x)) , (A.3) 

where = XB / '(fc^T) > is the inverse temperature if we choose the noise intensity 
D as in To check that flA.3h is the density of the invariant measure for (JT])-© and 
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([()]), one needs to verify that this p(x) is the solution of the stationary Fokker-Planck 
equation for ©-([9]), (Q. Such calculations are available, e.g. in [39|. 

Appendix B. 

In this Appendix we recall some generic facts from stochastic numerics [21]. In 
particular, we define the weak order of convergence of numerical methods for stochastic 
differential equations (SDEs) and discuss errors arising in computing ergodic limits. 
Let us introduce a system of SDEs of a general form 

r 

dX = a(X)dt + AP0dWi(*)> X(0) = x, (B.l) 
i=i 

where X, a, fii are <i-dimensional column- vectors and Wi(t), I = 1, . . . , r, are independent 
standard Wiener processes. Consider a numerical method for ( 1B.1I) based on the one- 
step approximation: 

X t:X (t + h)~ X ttX (t + h)=x + A(t, x,h;£), 0<t<t + /i<t*, (B.2) 

where £ is a random vector with moments of a sufficiently high order and A is a d- 
dimensional vector function. Introduce (for simplicity) the equidistant partition of the 
time interval [0,t*] into N parts with the step h = t±/N: = to < ti < • • • < 1>n — t*, 
£fc_i_i — t k = h. According to flB.21) . we construct the sequence 

X = x, X k+1 = X k + A(t k} X k , h; k = 0, . . . , N - 1, (B.3) 

where £i is independent of X and £, k+ i for k > is independent of X , . . . , X k , £i, . . . , 

We note that (1B.3j) contains both explicit and implicit one-step schemes. In 
explicit integration schemes the approximate solution at the next time-step, X k+ i, can 
be computed explicitly from the previous time-step value X k . For implicit methods, 
A(t, x, h; £) is a solution of an implicit relation with respect to x, i.e., implicit schemes 
in general require additional work. 

We usually distinguish two types of convergence of numerical methods for SDEs: 
mean-square (also called strong) and weak [24J. Mean-square methods are used for direct 
simulation of SDEs' trajectories which, e.g., can give information on general behavior 
of a stochastic model. Weak methods are sufficient for evaluation of mean values and 
are simpler than mean-square ones. We say that the method ( IB. 31) is weakly convergent 
with order p > if 

| (<p(X N )) - (<p(X(Q)) | < ChP (B.4) 

for functions (p which, together with their derivatives of a sufficiently high order, have 
growth at infinity not faster than polynomial. If a method converges with an order p in 
the mean-square sense, it also converges in the weak sense with order equal to or larger 
than p. The opposite is not true. Since weak methods suffice for computing averages, 
they are appropriate for the purposes of this paper. 
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To evaluate the expectation (<p(X N )) on a computer, one can apply the Monte 
Carlo technique: 

1 M 

m=l 

where X™ , m = 1, . . . , M, are independent realizations of the random variable X^- In 
(1B.5|) the first approximate equality involves the numerical integration error (cf. (1B.4j) ) 
and the error in the second approximate equality (the statistical error) comes from the 
Monte Carlo technique. 

The error of the Monte Carlo method in (1B.5I) is evaluated by 

[VarMXN)}] 1 / 2 

where, e.g., the values c = 1,2,3 correspond to the fiducial probabilities 0.68, 0.95, 
0.997, respectively, with the practical implication that 

u G (u -=\f%, u H ^=Vt)) , (B.6) 

1 M 



m=l 

with probability 0.68 for c = 1, 0.95 for c = 2, and 0.997 for c = 3. 

Now we assume that the solution of ( 1B.1I) is ergodic. In computing ergodic limits an 
additional error arises. We note that ergodic limits can be computed using the ensemble 
averaging or time averaging. In the former case it follows from a relation of the form 
(lA.ip for the solution X(t) of ( IB. II) that for any e > there exists t a > such that for 
all > t a 

\(<p(X x (Q)) - <p^\ < e. (B.7) 
Then we can use the following estimator for the ergodic limit (p er9 : 

AI 
m=l 

where the first approximate equality corresponds to the time cut-off while the second 
one relates to the numerical integration error, and the third to the statistical error as 
before. In this ensemble-averaging approach each of the errors is controlled by its own 
parameter (see [40J). 

The time-averaging approach to computing ergodic limits is based on a relation of 
the form ( )A.2j) . By approximating a single trajectory, one gets for a sufficiently large 

** L 
<f ~ i / <p(X x (s))ds ~ r 9 = jJ2^(X k ), (B.9) 

where Lh = Let us emphasize that £* in ( 1B.9j) is much larger than t* in (IB. 8|) because 
£* should be such that it not just ensures the distribution of X(t) to be close to the 
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invariant distribution (like it is required from £*) but it should also guarantee smallness 
of the variance of (p er9 . See further details about computing ergodic limits in, e.g. 
[40} UU 02] and the references therein. 

References 

Landau L D and Lifshitz E M. Phys. Zs. Sowjet, 8:153, 1935. 

Akhiezer A I, Bar'yakhtar V G, and Peletminskii S V. Spin Waves. North Holland, Amsterdam, 
1968. 

Vonsovsky S V. Magnetism. Wiley, New York, 1974. 

Aharoni A. Introduction to the Theory of Ferromagnetism. Oxford University Press, Oxford, 
2000. 

Zuti'c I, Fabian J, and Das Sarma S. Rev. Mod. Phys., 76:323, 2004. 

Gerrits Th, van den Berg HAM, Hohlfeld J, Bar L, and Rasing Th. Nature, 418:509, 2002. 
Kimel A V, Kirilyuk A, Usachev P A, Pisarev R V, Balbashov A M, and Rasing Th. Nature, 
435:655, 2005. 

Koopmans B, Ruigrok J J M, Dalla Longa F, and de Jonge W J M. Phys. Rev. Lett., 95:267207, 
2005. 

Melnikov A, Povolotskiy A, and Bovensiepen U. Phys. Rev. Lett., 100:247401, 2008. 
Novoselov K S, Gcim A K, Dubonos S V, Hill E W, and Grigorieva I V. Nature, 426:812, 2003. 
Dobrovitski V V, Katsnelson M I, and Harmon B N. Phys. Rev. Lett, 90:067201, 2003. 
Antropov V P, Katsnelson M I, van Schilfgaarde M, and Harmon B N. Phys. Rev. Lett., 75:729, 
1995. 

Antropov V P, Katsnelson M I, Harmon B N, van Schilfgaarde M., and Kuznecov D. Phys. Rev. 
B, 54:1019, 1996. 

Skubic B, Hellsvik J, Nordstrom L, and Eriksson O. J. Phys.: Condens. Matter, 20:315203, 2008. 
Hcllsvik J, Skubic B, Nordstrom L, Sanyal B, Eriksson O, Nordblad P, and Svedlindh P. Phys. 

Rev. B, 78:144419, 2008. 
Skubic B, Peil O E, Hellsvik J, Nordblad P, Nordstrom L, and Eriksson O. Phys. Rev. B, 79:024411, 
2009. 

Nowak U. In D. Stauffer, editor, Annual Reviews of Computational Physics IX, page 105. World 

Scientific, Singapore, 2001. 
Kazantseva N, Nowak U, Chantrell R W, Hohlfeld J, and Rebei A. Europhys. Lett., 81:27004, 
2008. 

d'Aquino M, Serpico C, Coppola G, Mayergoyz I D, and Bertotti G. J. Appl. Phys., 99:08B905, 
2006. 

Hairer E, Lubich C, and Wanner G. Geometric Numerical Integration: structure preserving 

algorithms for ordinary differential equations. Springer, 2002. 
Frank J, Huang W, and Leimkuhler B. J. Comp. Phys., 133:160, 1997. 
Arponen T and Leimkuhler B. BIT Num. Math., 44:403, 2004. 

Milstein G N, Repin Yu M, and Tretyakov M V. SI AM J. Numer. Anal, 40:1583, 2002. 
Milstein G N and Tretyakov M V. Stochastic Numerics for Mathematical Physics. Springer, 2004. 
Davidchack R L, Handel R, and Tretyakov M V. J. Chem. Phys., 130:234101-14, 2009. 
Krech M, Bunker A, and Landau D P. Comp. Phys. Comm., 111:1, 1998. 
Steinigeweg R and Schmidt H J. Comp. Phys. Comm., 174:853, 2006. 

Banas L. In Z. Li et al, editor, Lecture Notes in Computer Science: Numerical Analysis and its 

Applications, page 158. Springer, Berlin/Heidelberg, 2005. 
Kubo R and Hashitsume N. Prog. Theor. Phys. Suppi, 46:210, 1970. 
Brown W F. Phys. Rev., 130:1677, 1963. 

Gardiner C W. Handbook of Stochastic Methods. Springer, 2004. 



Stable and fast semi-implicit integration of the stochastic Landau- Lifshitz equation 24 

[32] Stratonovich R L. Conditional Markov Processes and their Application to the Theory of Optimal 
Control. Elsevier, 1968. 

[33] Hasminskii R Z. Stochastic Stability of Differential Equations. Sijthoff & Noordhoff, 1980. 

[34] Milstein G N and Tretyakov M V. J. Stat. Phys., 77:691, 1994. 

[35] Milstein G N and Tretyakov M V. SIAM J. Numer. Anal, 34:2142, 1997. 

[36] Shubin S and Zolotukhin M. Zh. Eksp. Teor. Fiz., 6:105, 1936. 

[37] Fisher M E. Am. J. Phys., 32:343, 1964. 

[38] Soize C. The Fokker- Planck Equation for Stochastic Dynamical Systems and its Explicit Steady 

State Solutions. World Scientific, 1994. 
[39] Garcia Palacios J L and Lazaro F J. Phys. Rev. B, 58:14937, 1998. 
[40] Milstein G N and Tretyakov M V. Physica D, 229:81, 2007. 

[41] Mattingly J C, Stuart A M, and Tretyakov M V. Convergence of numerical time-averaging and 

stationary measures via Poisson equations. arXiv:0908.4450v2, 2009. 
[42] Talay D. Stochastics and Stochastics Reports, 29:13, 1990. 



